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We simulate the sintering of particle aggregates due to surface diffusion. As a method we use Kinetic Monte-Carlo 
simulations in which elasticity can explicitly be taken into account. Therefore it is possible to investigate the 
shape relaxation of aggregates also under the influence of an external pressure. Without elasticity we investigate 
the relaxation time and surface evolution of sintering aggregates and compare the simulations with the classical 
Koch-Friedlander theory. Deviations from the theoretical predictions will be discussed. 



1 INTRODUCTION 

One interesting aspect of powder processing is ther- 
mal sintering: Due to atomic diffusion solid bridges be- 
tween particles form and grow. This lowers the surface 
free energy. If one waits long enough, complete coa- 
lescence of the particle aggregate will occur. The typi- 
cal time scales depend on the particle sizes and shapes, 
the material type and the prevailing temperature. Here 
we address open questions concerning the sintering dy- 
namics of large aggregates. 

It has long been known, that the time for the sin- 
tering of two identical particles is proportional to the 
4th power of the radius of the final particle. This 
is at least true for large particle s above the rough- 
ening temperature (Nichols 1966). For lower temper- 
atures the equilibra tion time increases exponentially 
JCombe et al. 20 00). 

In earlier investigations elasticity has not explicitly 
been taken into account, i.e. the atoms were restricted 
to discrete lattice positions. However, an elastic defor- 
mation changes the diffusion constant and hence in- 
fluences the sintering process. For example, the lat- 
tice constant of nano-particles may differ from its bulk 
value due to surface tension, or particles may be com- 
pressed differently in a powder under external load, de- 
pending on their position in the contact force network. 
Therefore we developed a program where the atoms 
are allowed to be displaced from their lattice sites in 
order to minimize the total elastic energy of the sys- 
tem. 

2 MODEL 

We use a three dimensional Kinetic Monte-Carlo 
(KMC) method to simulate the sintering of aggregates 



of nano-particles on a fee lattice. This means that grain 
boundaries are neglected in this paper. In our simula- 
tion the atomic displacements are calculated by finding 
the nearest local energy minimum by mea ns of a con- 
jugate gradient method (Pres s et al. 19921) . In order to 
save computing time, we update only the neighborhood 
of diffusing atoms, and relax the whole system elasti- 
cally only in the beginning and, if necessary, again af- 
ter rather large time intervals. The activation energies 
for the hopping rates then depend on the relaxed real 
positions of the atoms. 

The atoms interact via a Lennard- Jones potential. 
We calculate the binding energies E^i as the sum of 
the pair interactions up to the fourth nearest neighbor 
on the underlying fcc-lattice . 

In a diffusion step an atom hops from its initial po- 
sition (nominal lattice site i) to an unoccupied one, 
which is next to i on the underlying fcc-lattice. The 
activation energy for such a move is the difference be- 
tween the energy at the saddle point and the energy at 
the initial position of the atom. For the calculation of 
the saddle point energy E sp one has to distinguish two 
cases: The final site may be stable (at least three occu- 
pied neighbors) or unstable. In the latter case one can 
actually not find a local energy minimum at this site 
(see Fig[l]). Then the atom continues to move from this 
intermediate site to a randomly chosen stable final po- 
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Figure 1. 2d illustration for an atom that hops to a stable (a), respec- 
tively to an unstable intermediate position (b) from where it imme- 
diately continues to one of the neighboring stable sites (c). 
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Figure 2. Equilibration time for two particles of different size plotted 
versus the number of atoms (T = 800K). The radius of particle one 
is 10 atoms, the radius for particle 2 is varied from 4 to 20 as denoted 
in the bottom of the graph. ) 



sition next to it (with some probability it may actually 
jump back to its starting position i). 

Accordingly the saddle point energy for the first case 
is taken as 



E sp — E< 



sp,0 ■ 



(i) 



The first term is a constant parameter. The second one 
describes the strain dependence of the saddle point en- 
ergy. E£ i and E£ f are the strain dependent contribu- 
tions to the binding energy at the initial (T) and final 
('f ') atom positions. The average is taken to guarantee 
symmetry and a is an empirical parameter. This ansatz 
is justified by the finding that the binding as well as the 
saddle poi nt energy depend appro ximately linearly on 
the strain.( Schroeder & Wolf 1996 ) 

If the hop ends in an unstable, intermediate position, 
this is approximately regarded as the saddle point. We 
calculate the energy of the atom at the intermediate site 
giving it the average displacement of the stable atoms 
on the neighboring sites. This energy is taken as saddle 
point energy in the second case. 

In the KMC simulation every process is selected ac- 
cording to its corresponding rate which is calculated 
from the activation energy by: 



q = v exp(-/3E a ). 



(2) 



v is a fixed attempt frequency, (3 = with the 

Boltzmann factor and the absolute temperature T. 

Without elast icity E^^ corresponds to the bond 
counting model( Newman & Barkema 1999) and the 
strain dependent terms in ([!]) vanish. 

3 RESULTS 

3 . 1 Relaxation of two clusters 
Crystalline particles above their roughening tempera- 
ture are not faceted but round. If they have equal size, 
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Figure 3. The same data as in fig.|2|plotted versus the reduced radius. 

it is well known that the equilibration time, in which 
the two particles coalesce into a single one, is: 



:7V 4 / 3 



(3) 



where N oc r 3 is the total number of atoms. 

But what happens if the two diameters differ? In 
Fig.|2|the equilibration time is plotted versus the num- 
ber of atoms. If © were valid, the double logarithmic 
plot should give a straight line with slope 4/3. Instead 
we find that r ocr 4 with the reduced radius of the two 
particles, 
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(4) 



(see Fig. |3j . 



3.2 Surface evolution during sintering 
The surface evolution during the sintering process 
is usually described by th e Koch-Friedlander theory 
(Koch & Friedland erT990l : 



dA 
~dt 



A- A 



eq 



(5) 



Here A(t) is the surface of the agglomerate, A Qq the 
surface in equilibrium, r a relaxation time. However, 
this equation is interpreted and applied in several dif- 
ferent ways in the literature. Sometimes r and A Qq 
are viewed as constants, meaning the global relaxation 
time and the equilibrium surface area the whole aggre- 
gate will have after coalescence. In this case Eq. 
describes an exponential decay. 

In many cases r and A Qq are regarded as time de- 
pendent, local quantities. One assumes that the aggre- 
gate coarsens homogeneously by pairwise coalescence 
of particles. As explained above, the time constant r is 
then determined by the current radius of the constituent 
particles, which after n pairwise coalescence steps is 
2 n / 3 times the radius of the initial primary particles. 
The sintering becomes more slowly. 
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Figure 4. Snapshots during one simulation. The initial configuration 
consists of 12 clusters with 8372 atoms, T = 800K. 



Using a continuously varying 

r(t) ex (V/A(t)) 4 



(6) 



where V is the constant solid volume of the ag- 
gregate and determines A eq , the theory has recently 
been successfully applied to experimental sintering 
data of Ni-particles to calculate the activation energy 
of the relevant diffusion process, whereas the eval- 
uation with a constant r giv es unreasonable energy 
values. (Tsygan ovet al. 2004] ) 

What is missing, is a microscopic justification of 
Eq. 0. Therefore we simulated the sintering of ag- 
glomerates and measured the surface area. Typical 
snapshots are shown in Fig. |4| We fitted the surface 
area with the analytical solutions obtained for constant 
r or with Eq. (|6|), respectively. The constant A Qq was 
obtained as the average asymptotic surface area for 
large times. 

One can distinguish two different stages in which 
either a constant or a variable r gives a better descrip- 
tion of the surface dynamics. At early times when the 
particles successively merge together and the average 
radius increases, the fit with a variable r gives a bet- 
ter description (Fig.|3a)). When the average diffusion 
length does not increase significantly any more, a con- 
stant r is a better assumption (Fig. |3b)), but surpris- 
ingly its value is smaller than at the end of the early 
time regime with variable r. 

3.3 Influence of elasticity on sintering 
In order to understand the influence of an external 
stress we allowed elastic displacements of the atoms as 
described above. We set up a configuration of two par- 
ticles with periodic boundary conditions in x-direction. 
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Figure 5. Fit with solution of J5j for constant and increasing r at the 
beginning (a) and the end (b) of the sintering. 

The strain is controlled by the total system size in x- 
direction, which is kept constant during each simula- 
tion. 

In Fig.|6|an initial and final configuration are shown. 
One expects that the final configuration is reached 
faster under compression than under tension. We mea- 
sure the average squared radius along the x axis: 
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with 



1 N 1 N 

= ]vE^ and z ' = nT, 



(7) 



(8) 



The value R 2 (t) reflects how far away the system is 
from the thermodynamic equilibrium state. Its time de- 
pendence is plotted in Fig.Q It is scaled with the initial 
value, R 2 (t = 0), in order to eliminate the effect of the 
Poisson ratio for a better comparison of the curves for 
different initial strains. The initial configuration (upper 
part of Fig. |6l> are two spheres, whose surface is not 
in thermal equilibrium. Therefore initially the curves 
raise beyond 1 due to surface roughening. 

Looking at the relaxation behaviour we find that 
both systems evolve asymptotically similarly. How- 
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Figure 6. Sintering of two clusters with periodic boundary conditions 
in x-direction. An external stress (compressive/tensile) is applied by 
changing the system size in x-direction. 



ever, the relaxation process is faster for compressive 
strain which becomes obvious from the larger slope 
in the beginning of the coalescence. The equilibrium 
state is therefore reached earlier. The reason is that 
a > 1 in Eq. © for hopping diffusion in Lennard- 
Jones-Systems, and the binding of atoms on the sur- 
face be comes weaker under later al compression of the 
surface ( Schroeder & Wolf 1996 ). This implies that the 
activation energy decreases under lateral compression. 

The cylindrical configuration reached at the end of 
our simulation probably is only metastable: Due to the 
periodic boundary conditions it corresponds to an in- 
finitely long solid cylinder which should undergo a 
Rayleigh instability and split up into separate spheres. 
We think that thermal fluctuations in this direction are 
the reason why R increases again for larger times. 

4 CONCLUSIONS AND OUTLOOK 

We developed a 3d-KMC simulation program to model 
the sintering of nano-particles. The special feature of 
the program is that atoms are not fixed at their lattice 
sites in contrast to common KMC codes. This allows 
us to analyse the effect of elasticty on the sintering be- 
havior. 

The validation of the Koch-Friedlander theory leads 
to interesting results: The assumption of a constant 
r can only be applied if the characteristic diffusion 
length does not change. At the beginning of the sin- 
tering process of an agglomerate consisting of many 
clusters the characteristic length scale changes. In this 
case the assumption of a variable r gives a better de- 
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Figure 7. Average squared radius perpendicular to symmetry axis, 
scaled by initial value for the configuration shown in Fig.|6| 

scription of the surface evolution. 

We found that the equilibration time for two clus- 
ters of different size follows the r 4 power law with the 
reduced radius ©. 

The influence of elasticity for our set of parameters 
is only felt in the case of an external force. There it is 
found that compressive stress leads to faster relaxation. 
Without external stress it seems that the influence of 
the surface tension is negligible. This needs not always 
be the case, as our potential does not show strong sur- 
face tension. 

Therefore a next step would be to choose a potential 
that shows a stronger surface tension. So far grain 
boundaries were neglected as we used a continuous 
fee lattice. In the current development of the code 
these grain boundaries will be implemented. 
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